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ABSTRACT 

I investigate the use of Pearson's chi-square statistic, the Maximum Likeli- 
hood Ratio statistic for Poisson distributions, and the chi-square-gamma statistic 
(Mighell 1999, ApJ, 518, 380) for the determination of the goodness-of-fit between 
theoretical models and low-count Poisson-distributed data. I demonstrate that 
these statistics should not be used to determine the goodness-of-fit with data 
values of 10 or less. 

I modify the chi-square-gamma statistic for the purpose of improving its 
goodness-of-fit performance. I demonstrate that the modified chi-square-gamma 
statistic performs (nearly) like an ideal statistic for the determination of 
goodness-of-fit with low-count data. On average, for correct (true) models, the 
mean value of modified chi-square-gamma statistic is equal to the number of de- 
grees of freedom (z/) and its variance is 2z/ — like the distribution for v degrees 
of freedom. Probabilities for modified chi-square-gamma goodness-of-fit values 
can be calculated with the incomplete gamma function. 

I give a practical demonstration showing how the modified chi-square-gamma 
statistic can be used in experimental astrophysics by analyzing simulated X-ray 
observations of a weak point source (S/N^ 5.2; 40 photons spread over 317 
pixels) on a noisy background (0.06 photons per pixel). Accurate estimates (95% 
confidence intervals/limits) of the location and intensity of the X-ray point source 
are determined. 



^NOAO is operated by the Association of Universities for Research in Astronomy, Inc., under cooperative 
agreement with the National Science Foundation. 
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1. INTRODUCTION 



The goodness-of-fit between an observation of N data values, Xj, with errors, cxj, and a 
model, mj, can be determined by using the standard chi-square statistic: 

N 
i=l 

The standard chi-square statistic is the appropriate statistic to determine the goodness-of-fit 
whenever the errors can be described by a normal (a.k.a. Gaussian) distribution. 

Let us consider a more complicated situation where all the data values come from a 
pure counting experiment where each measurement^], n^, is a random integer deviate drawn 
from a Poisson ( |1837| , p. 205 et seq.) distribution, 

P{k■^,)^^e-\ (2) 

with a mean value of /i. The use of equation (|l]) in analyzing Poisson-distributed data is 
technically never correct. While the Poisson distribution approaches the normal distribution 
as the Poisson mean approaches infinity, a Poisson distribution never actually becomes a 
normal distribution even at very large Poisson mean values. The normal distribution is 
always symmetric; the coefficient of skewness for the normal distribution is zero. Poisson 
distributions are always asymmetric; the coefficient of skewness for a Poisson distribution of 
mean /i is While the Poisson distribution is almost symmetric about the mean for 

large mean values, its shape becomes progressively more asymmetric as the mean approaches 
zero. Thus the standard assumption that a Poisson distribution is approximately normally 
distributed is a good approximation only when the coefficient of skewness is negligible (e.g.. 

How does one then determine the goodness-of-fit with Poisson-distributed data? His- 
torically, many statistics have been proposed for the analysis of Poisson-distributed data. 
This paper will investigate the following four: 



a.: 



Pearson's x^'- 

where the expectation value of the mean of the parent Poisson distribution of the ith data 
value is assumed to be equal to the Poisson deviate [(/ij) = rii] and the square of the 



^ For example, X-ray photons, molecules, stars, galaxies, et cetera. 
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measurement error is assumed to be equal to the mean of the model Poisson distribution 



the modified Neyman's x ■ 



2 ^{rii-miY 

i= 



;|max(nj, Ij 

where the expectation value of the mean of the parent Poisson distribution of the ith data 
value is assumed to be equal to the Poisson deviate [( /Xi ) = nj\ and the square of the 
measurement error is assumed to be equal to Poisson deviate or one — whichever is greater 
[af = max(nj, 1)]; 



the Maximum Likelihood Ratio statistic for Poisson distributions: 

N 



i=l 



nii - rii + Uilia { — 



rrii 



(5) 



see, e.g., Baker & Cousins |1984] and references therein) 



and the chi-square-gamma statistic (Mighell |1999| ; hereafter PaperI): 



2 ^ + min {rii, 1) - rjiiY 

X-y — 2-^ 

i=l 



Tli + l 



(6) 



where the expectation value of the mean of the parent Poisson distribution of the ith data 
value is assumed to be equal to the Poisson deviate plus a small correction factor of zero 
for zero deviates and one in all other cases [(/Xj ) = rii + min(nj, 1)] and the square of the 
measurement error is assumed to be equal to the Poisson deviate plus one [af = + 1] . 

In PaperI, I demonstrated that the application of the standard weighted mean formula, 
Y.iiT'i'^r'^ I Hi^^i'^ 1 to determine the weighted mean of data, n^, drawn from a Poisson 
distribution, will, on average, underestimate the true mean by ~1 for all Poisson mean values 
larger than ~3 when the common assumption is made that the error of the ith observation 
is (Tj = max(y^, 1). This small, but statistically significant offset, explains the long-known 
observation that chi-square minimization techniques which use the modified Neyman's 
statistic [eq. (^)] to compare Poisson-distributed data with model values, mj, will typically 
predict a total number of counts that underestimates the actual total by about 1 count per 
bin (see, e.g., Bevington |196|, Wheaton et al. |T995|) . 



Based on my finding that the weighted mean of data drawn from a Poisson distribu- 
tion can be determined using the formula X^i ['^j + niin ip-i-, 1)] (^j + / Hi {^i + 1)~ 
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I proposed that the chi-square-gamma statistic, [eq- (0)], should always be used to analyze 
Poisson-distributed data in preference to the modified Neyman's statistic. Following my 
own advice, I will not discuss the modified Neyman's statistic in the remainder of this 
article. 

The chi-square distribution for z/ degrees of freedom approaches a Gaussian distribution 
with a mean equal to u (i.e., ^ = v) and a variance equal to 2v (i.e., = 2v) as the number 
of degrees of freedom approaches infinity. Ideally, a x^ statistic for Poisson distributions 
for V (independent) degrees of freedom would exhibit the same behavior as the number of 
degrees of freedom approaches infinity for all Poisson mean values (i.e., /i>0). 

Do the Xp; x\i a.nd the x^ statistics perform as expected for large Poisson mean val- 
ues? These three x^ statistics are applied to the same data set in Fig. |l| [top to bottom, •fig0 
respectively). For this example, an ideal x^ statistic for Poisson-distributed data would 
have a cumulative distribution similar to that of the chi-square distribution for 10"^ degrees 
of freedom which well approximated as the cumulative distribution function of a Gaussian 
distribution with a mean of 10^ and a variance of 2x 10"^. The results of the top and bottom 
panels are well matched to the expected cumulative distribution; the differences between the 
expected and measured mean and rms values are not statistically significant. The Xp and 
the X7 statistics perform as expected with a Poisson mean value of 100. The cumulative 
distribution of the middle panel, however, clearly deviates from the expected cumulative 
distribution; the difference between the expected and measured mean and rms values, while 
small, is statistically significant. The Xa does not perform like an ideal x^ statistic for Pois- 
son distributions with a mean value of 100 — a level that is generally considered to be well 
above the low-count regime (/i <^ 25). 

Let us continue the investigation of the performance of the Maximum Likelihood Ratio 
statistic for Poisson distributions with 1000 samples of 10^ Poisson deviates with Poisson 
mean values of 100, 10, 1, 0.1, and 0.001. Figure || confirms that the x\ statistic does not •fig| 
perform like an ideal x^ statistic in the low-count regime. The average contribution by the 
ith deviate to an ideal x^ statistic for the analysis of Poisson-distributed data would be 
exactly one and the average contribution to its variance would be exactly two. Figure ^ •fig| 
expands the previous analysis of PaperI of the Xa statistic over a wide range of Poisson 
mean values from 0.001 to 1000. The dashed lines of Fig. ^ show the results for an ideal x^ 
statistic; one can clearly see that the average contribution to Xa is not equal to one and the 
average contribution to its variance is not equal to two for Poisson mean values <10. The 
poor performance of the Xa statistic with low-count data may come as a surprise to many 
readers since it has historically been advocated as being one of the best statistics for the 
analysis of Poisson-distributed data. 



- 5 - 



Chi-square statistics can serve (at least) two distinct purposes: (1) their functional 
forms can be utilized as the core of parameter estimation algorithms, and (2) their values 
can serve as a measure of the goodness-of-fit between a model and a data set. 
• While the functional form of the Xx statistic can be successfully utilized for the purpose 
of parameter estimation with Poisson-distributed data in the low-count regime (see, e.g., 
PaperI), the Maximum Likelihood Ratio statistic for Poisson distributions [eq. ^] should 
not be used to determine the goodness-of-fit with low-count data where the Poisson mean is 

In this work, I investigate the use of Pearson's statistic and the chi-square-gamma 
statistic for the determination of the goodness-of-fit between theoretical models and data 
derived from counting experiments. I develop a methodology in §2 which modifies Pear- 
sons's chi-square statistic for the purpose of improving its goodness-of-fit performance. This 
methodology is then be applied to modify the chi-square-gamma statistic (§3). The modi- 
fied chi-square-gamma statistic is shown to perform (nearly) like an ideal statistic for the 
determination of goodness-of-fit with low-count data. Simulated X-ray images are analyzed 
in §4 as a practical demonstration of the possible use of the modified chi-square-gamma 
statistic in experimental astrophysics. The summary of the paper is presented in §5. 



2. THE MODIFIED PEARSON'S STATISTIC 

Let us continue the investigation of the performance of Pearson's x^ statistic with 1000 
samples of 10^ Poisson deviates with Poisson mean values of 100, 10, 1, 0.1, and 0.001. 
Figure ^ shows that the Xp statistic does not perform like an ideal x^ statistic for Poisson •fig| 
mean values <^ 10. The average contribution by the ith deviate to an ideal x^ statistic for 
the analysis of Poisson-distributed data would be exactly one and the average contribution 
to its variance would be exactly two. Figure ^ expands the previous analysis of PaperI of •fig| 
Pearson's x^ statistic over a wide range of Poisson mean values from 0.001 to 1000. The 
dashed lines of Fig. || show the results for an ideal x^ statistic; one can see that while the 
average contribution to Xp is one, the average contribution to its variance is not equal to 
two for Poisson mean values ^10. 

• Pearson's x^ statistic [eq. ^] should not be used to determine the goodness-of-fit with 
low-count data where the mean of the parent Poisson distribution is < 10. 

The variance of Pearson's x^ statistic is, by definition, 

N 
i=l 



j=i 'J 



(7) 
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where u = N — M is the number of independent degrees of freedom, N is the number of data 
values, and M is the number of free parameters. The variance of the reduced chi-square of a 
statistic for a large number of observations should ideally be two. In the limit of a large 
number of observations of a single Poisson distribution with a mean value of fi, the variance 
of the reduced chi-square of the Pearson's statistic is 



lim 



{rii - rriif 



{Ui - ruif 



{Ui - niif 




1 ^ {rij-mjf 









n 


V 





lim 



Ui — rriA 



rrii 



Hp 



- 1 



Ui - lim [fip]) 

N^oo J ^ 



lim [^p 



rij - hY 



^ lim 



N 



-J2{NPik;fi)} 

r ; n 



k=0 



1 oo 

k=Q 
1 

2 + - . 



[see eq. (25) of PaperI] 
[see eq. (5) of PaperI] 



[see eq. (7) of PaperI] 



(8) 



If we assume that Pearson's applied to a large number of observations of a single 
Poisson distribution with a mean value of /z always produces a normal distribution with a 
mean equal to the number of degrees-of-freedom (z/) [see eq. (25) of PaperI] and a variance 
of z/(2 + /i~^) [see eq. (P)], we can then attempt to create an ideal statistic for the analysis 
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of Poisson-distributed data by modifying Pearson's as follows: 



N 



XpM — X! 



i=l 



{xli - ixh)} 



1/2 



+ 1 



where 



xli = 



I Hi - nii 



rrii 



(9) 



(10) 



is the contribution of the ith data value to Pearson's x^? 



ix 



2 ] 
Pi) 



IS 



the expectation value of Xpi [see eq. (25) of PaperI], 



2 + m- 



(11) 



(12) 



is the variance of (xpj) [see eq. (H)]. Translating the mathematical notation to English, we 
have (1) shifted the mean of the standard Xp distribution from u times equation (|1T]) to zero, 
(2) forced the variance of the shifted distribution to be exactly 2z/, and then (3) shifted the 
mean of the variance-corrected distribution from zero back to u. Thus, by definition, the 
modified Pearson's chi-square statistic (xpm) ^ill have a mean value of u and a variance of 
2z/ — in the limit of a large number of observations. 

Let us now investigate the performance of the modified Pearson's statistic with 1000 
samples of 10^ Poisson deviates with Poisson mean values of 100, 10, 1, 0.1, and 0.001. 
Figure || shows that Xpm results are significantly better than Xp results [Fig. ^ — especially 
for Poisson mean values less than 10. Figure |^ investigates the performance of the modified 
Pearson's x^ statistic over a wide range of Poisson mean values from 0.001 to 1000. The 
dashed lines of Fig. ^ show the results for an ideal x^ statistic; one can see that while the 
average contribution to Xpm is 1, as expected, and the average contribution to its variance 
is equal to 2, as expected, the performance is not uniform for all Poisson mean values — the 
spread seen in the variance plot {bottom panel) increases as the Poisson mean approaches 
zero. 

Figures |^ and indicate the the modified Pearson's x^ statistic works well in the perfect 
case where one has a priori knowledge of the true Poisson mean. In a real experiment, the 
true mean of the parent Poisson distribution is rarely (if ever) known and model parameters 
must be estimated from the observations. How well does the modified Pearson's x^ statistic 
work with reasonable parameter estimates? Comparing Fig. |^ with Fig. ^ and Fig. ^ with 
Fig. 0, we see that the variances are significantly smaller when a realistic model (i.e., the 
sample mean) is used instead of a perfect model (i.e., the true mean). A statistic that fails 



•Fig| 



• FlG^ 
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with reasonable parameter estimates is not a very useful statistic for the analysis of real 
observations with low-count data. 

• The modified Pearson's statistic [eq. ij^] should not be used to determine the goodness- 
of-fit with low-count data where the mean of the parent Poisson distribution is ^10. 



3. THE MODIFIED CHI-SQUARE-GAMMA STATISTIC 

Let us continue the investigation of the performance of the chi-square-gamma statistic 
with 1000 samples of 10^ Poisson deviates with Poisson mean values of 100, 10, 1, 0.1, and 
0.001. Figure |TU| confirms that the statistic does not perform like an ideal statistic in 
the low-count regime. Figure expands the previous analysis of PaperI of the statistic 
over a wide range of Poisson mean values from 0.001 to 1000. The X7 statistic clearly does 
not perform like an ideal statistic for Poisson mean values < 10. 

• The chi-square-gamma statistic [eq. (Qj/ should not be used to determine the goodness- of- fit 
with low-count data where the mean of the parent Poisson distribution is ^10. 



•Fig|10| 
•Fig|ii| 



The variance of the chi-square-gamma statistic is, by definition. 



0-^2 



N 

E 



[rii + min(nj, 1) 



mi 



rii + l 



V 



N 

E 



[Uj + min(nj, Ij 



rij + 1 



(13) 

-i=i -J^' 

where v = N — M is the number of independent degrees of freedom, N is the number of data 
values, and M is the number of free parameters. The variance of the reduced chi-square of a 
statistic for a large number of observations should ideally be two. In the limit of a large 
number of observations of a single Poisson distribution with a mean value of ^, the variance 
of the reduced chi-square of the chi-square-gamma statistic is 
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Ui + l 

[ui + min(nj, 1) — rrii 



V 



rij + 1 
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E 



[ui + min(?T,j, 1) — rrii 
rii + l 



[see eq. (29) of PaperI] 
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lim 



lim 

N—>oo 



lim 

N-ioo 



lim 

N^oo 



N 



N 



1^ 

^ 1=1 



[rii + min(nj, 1) — yU^] 
rii + l 



1 + e-^" ifi - 1) 



[see eq. (18) of PaperI] 



TV 



N 



-\ 2 



Ui + min(r;,j, 1) — lim [/x^] 



i=l 



N 



^ i=l 



rii + l 



[rii + min(rai, 1) - 



1 + e-^ (/i - 1) 



1 



Ui + l 

'' [k + min(fc, 1) — fif 



1 + e-'^ (/i - 1) 



[see eq. (19) of PaperI] 



fc=0 



fc=0 



[k + min(fc, 1) — fiY 
k + 1 



k + 1 



1 + e-^' (/i - 1) 



l + e^^ (/X- 1) 



/i^e"'' [Ei(/i) - 7em - ln(/i) + 4] - /i^ - /i + e"^ -2/2^ + 2/i + 1 



+ e 



-2^ 



+ 2/i - 1 , (14) 



where Ei(x) is the exponential integral of x [Ei(x) 
and 7em = lim^^oo {j27=i ^} - ln(n) 



oo t 



■ dt for X > 0] 



i 0.5772156649 is the Euler-Mascheroni constant. 
Equation ([T^) approaches the expected value of 2 for large Poisson mean values [see the 
solid curve in the bottom panel of Fig. Oj . 



If we assume that chi-square-gamma statistic is applied to a large number of observations 
of a single Poisson distribution with a mean value of always produces a normal distribution 
with a mean equal to the number of degrees of freedom (z/) times see equation (29) of PaperI 
and a variance of v times equation ([T^) , we can then attempt to create an ideal statistic 
for the analysis of Poisson-distributed data by modifying the chi-square-gamma statistic as 
follows: 



X7M 



N 

E 

i=l 



1/2 



where 



2 + ™^ i^h 1) - 



rii + l 

is the contribution of the ith data value to the chi-square gamma statistic, 



l + e 



[mi 



1) 



(15) 

(16) 
(17) 



is the expectation value of x^i [see equation (29) of PaperI], and 



a 



m. 



|e [Ei(mi) - 7em - In(mi) + 4] - - + 



-2ml + 2m j + 1 



+ e 



-2mi 
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is the variance of (x^j) [see equation (n)]. Translating the mathematical notation to English, 
we have (1) shifted the mean of the standard Xp distribution from z/ times equation (|1^) to 
zero, (2) forced the variance of the shifted distribution to be exactly 2z/, and then (3) shifted 
the mean of the variance-corrected distribution from zero back to u. Thus, by definition, the 
modified chi-square statistic statistic (X7M) ^^^^ have a mean value of u and a variance of 2z/ 
— in the limit of a large number of observations. 

Let us now investigate the performance of the modified chi-square-gamma statistic with 
1000 samples of 10^ Poisson deviates with Poisson mean values of 100, 10, 1, 0.1, and 0.001. 
Figure |12| shows that X7M results are significantly better than the X7 results [Fig. |T0[ — 
especially for Poisson mean values less than 10. Figure O investigates the performance of 
the 

X^M statistic over a wide range of Poisson mean values from 0.001 to 1000. The dashed 
lines of Fig. O show the results for an ideal statistic; one can see that while the average 
contribution to x'^u is 1, as expected, and the average contribution to its variance is equal 
to 2, as expected, the performance is not uniform for all Poisson mean values. The bump 



seen in the bottom panel of Fig. ^ near the Poisson mean value of 10 is an artifact caused 
by the min(ni, 1) offset in the numerator of the definition of the chi-square-gamma statistic 
[eq. 



Figures |I2| and |I3] indicate the the modified chi-square-gamma statistic works well in the 
perfect case where one has a priori knowledge of the true Poisson mean. How well does the 
modified chi-square-gamma statistic work with reasonable parameter estimates? Comparing 
Fig. with Fig. ^ and Fig. |T5| with Fig. |l^, we see that the results for the modified X7 
statistic with a realistic model (i.e., the sample mean) are nearly identicalQto those obtained 
with a perfect model (i.e., the true mean). 

• The modified chi-square-gamma [eq. (|7^/ statistic performs (nearly) like an ideal x^ 
statistic for the determination of the goodness- of- fit with low-count data. On average, for 
a large number of observations, the mean value of x'^u statistic is equal to the number of 
degrees of freedom (u) and its variance is 2v — like the x^ distribution for v degrees of 
freedom. 



•Fig|12| 

•Fidii 



•Fig 
•Fig 



14 
15 



^ The measured mean values appearing on the right side of top 3 panels of Fig. 14 are about 1 lower 
than the comparable value given in Fig. |l2|. This is the result of losing one degree-of-freedom due to the 
determination of the sample mean from the data (i.e., v drops from 10000 to 9999). The scrambling caused 
by the modification of the statistic appears to have caused this expected loss of one degree-of-freedom to 
vanish in the very-low-count data regime (/i^O.l). 
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4. SIMULATED X-RAY IMAGES 

I now demonstrate the new modified chi-square-gamma statistic by using it to study 
simulated X-ray images. Cash ( |1979| ) apphed his C statistic to the problem of determining 
the position of a weak source in a X-ray image. Let us use Cash's Point Spread Function 
(PSF) but with a resolution of 100 pixels per unit area: 

{vr 
— (1 — r) for r < 1, , ^ 

3^ \ - ' (19) 
for r > 1, 

where = (x/10)^ + (y/10)^. This PSF has the volume integral of 

„3' 




' ^ (20) 
for r > 1. 



Figure |T6| shows a simulated observation of a point source with an intensity of 40 X-ray •fig|i6| 
photons on a background flux of 0.06 X-ray photons per pixel. This observation contains 
2786 pixels with photons, 204 pixels with 1 photons, and 10 pixels with 2 photons. There 
are a total of 56 photons found in the 317 pixels within a radius of 10 pixels of the center 
of the X-ray point source which is located at the (x, y) position of (33, 26). This is clearly a 
marginal detection of a weak X-ray point source on a noisy background; the peak signal-to- 
noise ratio (~5.2) occurs at a radius of ~8 pixels. 

We will now use the modified chi-square-gamma statistic to answer the following ques- 
tions about this X-ray image: 

1. Is there an X-ray point source in the image? 

2. If so, where is it located? 

3. What is its total intensity? 

The exact determination of the location and intensity of the X-ray point source in Fig. ITB 



is precluded by the fact that this particular observation contains only low-count data — we 
must be content with realistic estimates for the location and intensity based on a detailed 
statistical analysis of the data. 

Our first objective is to determine if there is an X-ray point source in the observation. 
One way this can be done is to investigate the region(s) containing non-point- source pixels 
(data values). This approach requires knowledge of the background flux level — which we 
will henceforth assume is constant throughout the entire image. We begin by making a 
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rough first estimate of tlie background flux level by dividing the total number of photon in 
the image by the total number of pixels: 0.0747 (~ 224 / 3000) photons per pixel. 

The background flux level estimate may be significantly improved with a bit more work. 
There are 18 photons in the 317 pixels within 10 pixels of the position (12, 15) of Fig. |16|. Is 
the detection of 18 photons consistent with the expected value of 23.6799 (= 0.0747 x 317) 
photons? The upper and lower 99.9% single-sided confidence limits for 18 photons are 35.35 
and 7.662, respectively [see Tables 1 and 2 of Gehrels |1986|| . I conclude that the pixel at 



(12, 15) is a background pixel (shown as such with a gray box in Fig. because the expected •fig|i^ 
number of background photons lies within the range of the upper and lower 99.9% single- 
sided confidence limits of the observed mimbei of photons (i.e., 7.662 < 23.6799 < 35.35). 
There are 56 photons in the 317 pixels within 10 pixels of the position (33, 26) of Fig. |IB[ 
The upper and lower 99.9% single-sided confidence limits for 56 photons are 83.1784 and 



35.6834, respectively [see eqs. (10) and (14) of Gehrels |1986|| . I conclude that the pixel at 
(33, 26) is not a background pixel because the expected number of background photons is 
less than the lower 99.9% single-sided confidence limit of the observed number of photons 
(i.e., 23.6799 < 35.6834). This conclusion was expected since (33,26) is the center of the 
X-ray source. The 2676 gray pixels in Fig. |1^ have a total of 162 photons. We can now 
make a second estimate of the background flux: 0.0605 (~ 162 / 2676) photons per pixel. 
Repeating this process once more yields the final estimate of the background flux: 0.0601 
(^ 153 / 2547) photons per pixel. The measurement error for this estimate is approximately 
0.0049 (^ V153 + 1 / 2547). The final estimate of the X-ray background flux, 0.0601±0.0049, 
is in excellent agreement with the true value of 0.06 [see Fig. •fig|i8| 



I conclude that Fig. has at least one X-ray point source because the entire data set 
is not consistent with a X-ray background flux of 0.0601 photons per pixel for every pixel. 
Assuming that there is only one X-ray source, we can make the first rough estimate of its 



location by stating that it probably is located at a non-gray pixel location in Fig. Il8l There 



are 453 non-background pixels in Fig. ITB| with a total of 71 photons. This fact allows us to 



restrict the uncertainty of the location of the X-ray point source to about 15% (~ 453 / 3000) 
of the total image. Assuming a background flux of 0.06 photons per pixel, we expect that 27 
{k, 453 X 0.06) photons of the total 71 photons found in the non-background pixels would be 
due to the background and not the point source. We can now make the first rough estimate 
of the intensity of the X-ray source: 44 (=71 — 27) photons. 

There are 48 photons in the 317 pixels within 10 pixels of the position (24,26) of Fig. 
IB]. Is this photon sum consistent with a model of a 40 photon point source centered at 



that location on a background of 0.06 photons per pixel? The upper and lower 95% single- 
sided confidence limits for 48 photons is 61.05 and 37.20 photons, respectively. The model 
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predicts that we should find 58.9802Q photons within a radius of 10 pixels. The 48 photons 
found within 10 pixels of (24, 26) is consistent with the model (shown as such with a dark- 
gray circle in Fig. |T^), because the expected number of photons lies within the range of 
the upper and lower 95% single-sided confidence limits of the observed number of photons 
(i.e., 37.20 < 58.9802 < 61.05). There are 217 circled pixels in Fig. [l|. This fact allows 
us to further restrict the uncertainty of the location of the X-ray point source to ~7.2% 
(fti 217 / 3000) of the total image. 

Since the peak signal-to-noise ratio occurs near a radius of 8 pixels, we now investigate if 
we can improve our estimate of the location of the X-ray point source by considering photon 
sums within a smaller aperture with a radius of 8 instead of 10 pixels. There are 34 photons 
in the 197 pixels within a radius of 8 pixels of the position (26, 26) of Fig. |TB|. Is this photon 



sum consistent with a model of a 40 photon point source centered at that location with a 
background of 0.06 photons per pixel? The upper and lower 95% single-sided confidence 
limits for 34 photons is 45.27 and 25.01 photons, respectively. The model predicts that we 
should find 47.2664 photons within a radius of 8 pixels. The 34 photons found within 8 pixels 
of (26, 26) is not consistent with the model since the expected number of photons is greater 
than the upper 95% single-sided confidence limit of the observed number of photons (i.e., 
47.2664 > 45.27). However, the 39 photons found within 8 pixels of (27,26) is consistent 
with a 40 photon point source centered at (27, 26) on a background of 0.06 photons per pixel 
(i.e., 29.33 < 47.2664 < 50.94). There are a total of 111 circled pixels in Fig. |T| . This .fig|i| 
fact allows us to further restrict the uncertainty of the location of the X-ray point source to 
~3.7% (= 111/3000) of the total image. 

One way to boost the data out of the low-count regime is to compare the cumulative 
radial distribution of the model with the cumulative radial distribution of the data. The 
modified chi-square-gamma statistic was used to compare the cumulative radial distribution 
of the model (10 1-pixel-wide bins =^ 10 degrees-of-freedom) with the cumulative radial 
distribution of the data (similarly formatted). At the position of (27,26) the value of x^m 
for the cumulative radial distributions was computed to be 13.7338 (z/ = 10) for a model 
of a 40 photon point source centered at (27, 26) on a background of 0.06 photons per pixel. 
The 95th percentage point for the chi-square distribution with 10 degrees of freedom may 
be found in several standard references: 18.3 {CRC Handbook of Chemistry and Physics, 



These computations only included pixels within an aperture if the center of the pixel was within the 
given aperture radius; partial pixels whose center was just outside aperture boundary were rejected. The 
minimum number we would expect the model to predict is 58.8496 [40 + (tt x 10^ x 0.0600)] photons. The 
maximum number we would expect the model to predict is 59.0200 [40+ (317 x 0.0600)] photons. The model 
prediction lies within these extremes: 58.8496 < 58.9802 < 59.0200. 
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Lide & Frederikse |T995|, p. A-106), 18.31 (Bevington |T96|, p. 315) and 18.3070 (Abramowitz 



& Stegun |1964| , p. 985). I conclude that the image location (27,26) is within the 95% 
confidence interval because the value of the modified chi-square-gamma statistic is less than 
the 95th percentage point for the chi-square distribution with 10 degrees of freedom (i.e. 
13.7338 < 18.3070). The probability that the observed chi-square value for a correct model 
should be less than a value of for ^ degrees of freedom is P(x^|i^) = -P(|, ^) where the 
latter function is the incomplete gamma function [ P = 1 — Q; see, e.g., the GAMMP routine 
in Numerical Recipes (Press et al. |1986|) j. If we assume that x^u is distributed like the x^ 
distribution, then we can assign a probability for the modified chi-square-gamma value for 
10 degrees of freedom: P(13.7338|10) = P(f , i^^) = 0.814517. There is thus a -81.5% 
chance that the observed modified chi-square-gamma statistic will be less than 13.7338 for 
10 degrees of freedom. The contour in Fig. |19| shows the 95% confidence interval of the 
X-ray point source based on the x^m analysis of the cumulative radial distribution of the 
data. The value of x^u fo^ the cumulative radial distribution at (26, 26) was computed to 
be 30.4707 giving a probability of ~99.9%; this location in Fig. [T^ hes outside the 95% 
confidence interval. 

We can further use the modified chi-square-gamma statistic with the cumulative radial 
distribution to determine the 95% confidence limits of the intensity of the X-ray source 
in the image (see Fig. EH ). The upper and lower single-sided 95% confidence limits for •figEoI 



the intensity of an X-ray point source at (33, 26) in Fig. |T6| is 54.5 and 28.0, respectively. 
The true intensity of the X-ray source is 40 photons. Given a background fiux uncertainty 
of as = 0.0049 photons per pixel (see above), we can approximate the theoretical rms 
measurement error for a 40 photon point source spread over 317 pixels {A = 317 px^) as 
a ~ V40 + 1 + A(Tb ~ 8.0 photons. The difference between the upper and lower 95% single- 
sided confidence limits is approximately 3.3 standard deviations of the normal probability 
function. This fact can be used to approximate an rms measurement error for our intensity 
estimate of ai ^ 8.0 (54.5 - 28.0)/(2 x 1.65)] photons. The 

X7M analysis using the 
cumulative radial distribution has yielded an excellent intensity estimate. 



The analysis presented in Figures and [2^ is predicated on the assumption that the 
modified chi-square-gamma statistic is distributed like x^- But is this assumption valid? 



Figure 21 shows that the analysis of 10^ simulated X-ray observations like Fig. yields •fig|i| 
modified chi-square-gamma values that are distributed like the chi-square distribution for 
z/ = 317 degrees of freedom: a Gaussian distribution with a mean of u and a variance of 
2z/. The above analysis has assumed that the probability that the observed modified chi- 
square-gamma value for a correct model should be less than a value of for degrees of 
freedom can be given as P{x'^mW)- Assuming that the predicted probability from P(x^mI^) 
is an accurate prediction of the true probability, then the predicted probability of the 9500th 



- 15 - 



simulated observation of a total of 10000 (sorted by x^m value) should be very close to 95%; 
Fig. indicates that this is indeed the case (i.e., the probability for the x^m value of the 
9500th simulated observation is 94.8666%). Since x^m statistic is distributed (nearly) like 
the x^ distribution, the usage of the incomplete gamma function to predict probabilities for 
modified chi-square-gamma values appears to be justified in practical analysis problems. 

5. SUMMARY 

I investigated the use of Pearson's chi-square statistic [eq. (|^)], the Maximum Likeli- 
hood Ratio statistic for Poisson distributions [eq. (^], and the chi-square-gamma statistic 
[eq. d^)] for the determination of the goodness-of-fit between theoretical models and low- 
count Poisson-distributed data. I concluded that none of these statistics should be used to 
determine the goodness-of-fit with data values of 10 or less. 

I modified Pearson's chi-square statistic for the purpose of improving its goodness-of-fit 
performance. I demonstrated that modified Pearson's x^ statistic [eq. (^] works well in 
the perfect case where one has a priori knowledge of the correct (true) model. In a real 
experiment, however, the true mean of the parent Poisson distribution is rarely (if ever) 
known and model parameters must be estimated from the observations. I demonstrated 
that the modified Pearson's x^ statistic has a variance that is significantly smaller than that 
of the x^ distribution when realistic models, defined as having parameters estimated from 
the observational data, are compared with Poisson-distributed data. Any statistic that fails 
with models based on reasonable parameter estimates is not a very practical statistic for the 
analysis of astrophysical observations. I concluded that the modified Pearson's x^ statistic 
should not be used to determine the goodness-of-fit with low-count data values of 10 or less. 

I modified the chi-square-gamma statistic for the purpose of improving its goodness- 
of-fit performance. I demonstrated that the modified chi-square-gamma statistic [eq. ([T5|)] 
performs (nearly) like an ideal x^ statistic for the determination of goodness-of-fit with 
low-count data. On average, for correct (true) models, the mean value of the modified chi- 
square-gamma statistic is equal to the number of degrees of freedom (z/) and its variance is 
2u — like the x^ distribution for v degrees of freedom. 

An ideal x^ statistic for the determination of goodness-of-fit with low-count data should 
fail in a predictable manner. Hypothesis testing of low-count Poisson-distributed data with 
the modified Pearson's x^ statistic will produce the peculiar and undesirable result that 
correct models are more likely to be rejected than realistic models [cf. Fig. |^ with Fig. ||. 
The modified chi-square-gamma statistic is a practical statistic to use for hypothesis testing 
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of astrophysical data from counting experiments because it performs (nearly) like an ideal 
statistic for realistic and correct models in the low-count and the high-count data regimes; 
accurate and believable probabilities for x^m goodness-of-fit values can be calculated with 
the incomplete gamma function [Figs. ^ and A lot of nothing can tell you something 
— as long as there are some observations with signal in them. 

Vincent Eke sent me an e-mail asking if I had an expression for the variance of the 
statistic which described the mysterious second hump of the solid curve of Fig. 3 of PaperI. 
After a rapid exchange of email with me over the period of a week, he was the first to derive 
an analytical formula for cr^2 [eq. ([l^)]- The knowledge that the variance of could in 
fact be expressed explicitly as an analytical expression turned out to be the breakthrough 
that 1 had needed in order to complete the development of the modified statistic. It is a 
pleasure to acknowledge his contribution to this research. 

I would like to thank Mike Merrill for the use of his copy of Mathematica which I used 
to check some of the arithmetic of the critical last step in the derivation of Eq. (0). 
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Fig. 1. — A simulated data set of 1000 samples ("observations") of 10"^ Poisson deviates 
("measurements") per sample was created assuming a mean value /i=100 for each Poisson 
deviate. Each sample in this data set was then analyzed using Pearson's statistic [top; 
definition: eq. (^)], the Maximum Likelihood Ratio statistic for Poisson distributions [mid- 
dle; definition: eq. (H)], and the chi-square-gamma statistic [bottom; definition: eq. (^]. The 
model of the ith deviate in each sample was set to the true mean value of parent Poisson 
distribution (i.e., rrii = fi = 100) and the number of independent degrees-of-freedom was 
therefore equal to the number of deviates per sample (i.e. u = 10^). Compare the cumula- 
tive distribution for each statistic with the cumulative distribution function of a Gaussian 
distribution with a mean of 10^ and a variance of 2 x 10^ [thick curve in each panel]. The 
number and error shown on the right side of each panel is the mean and rms value of the 
1000 samples shown in that panel; ideally these values should be about 10000±141.4. 
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Fig. 2. — The cumulative distribution functions for 1000 samples of 10^ Poisson deviates 
{top to bottom: fi = 100, 10, 1, 0.1, and 0.01) analyzed using the Maximum Likelihood Ratio 
statistic for Poisson distributions [definition: eq. (|^)]. In all cases, i/ = 10^ and was set to 
the true mean value of the data set. Other details as in Fig. |1|. 
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Fig. 3. — Top panel: Reduced chi-square as a function of the true Poisson mean (0.001 < 
fi < 1000 with 10 mean values per decade) for the Maximum Likehhood Ratio statistic 
for Poisson distributions with the model of the ith deviate set to the mean value of the 
parent Poisson distribution. The open squares show the results of the analysis of one sample 
composed of 10^ Poisson deviates (z/ = 10^) at each given Poisson mean value. The filled 
squares show the results of the analysis of 1000 subsamples of the lO'' Poisson deviates 
(z/ = 10^) previously analyzed as one large sample. The scatter of the filled squares with 
respect to the open squares is real and is due to random fluctuations of the parent Poisson 
distributions. The dashed line shows the ideal value of one. Bottom panel: The variance of 
the reduced chi-square values shown in the top panel. The dashed line shows the ideal value 
of two. 
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Fig. 4. — The cumulative distribution functions for 1000 samples of 10^ Poisson deviates 
[to]) to bottom: fi = 100, 10, 1, 0.1, and 0.01) analyzed using Pearson's statistic [definition: 
eq. @] (same input data set as for Fig. H). In all cases, z/ = 10^ and was set to the true 
mean value of the data set. Other details as in Fig. 




Fig. 5. — Reduced chi-square as a function of the true Poisson mean for Pearson's 
statistic with the model of the ith deviate set to the true mean value of the parent Poisson 
distribution (same input data set as for Fig. |^). The solid line connecting the open squares 
in the bottom panel is the formula 2 + [see eq. (^]. Other details as in Fig. ^. 
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Fig. 6. — The cumulative distribution functions for 1000 samples of 10^ Poisson deviates 
{top to bottom: fi = 100,10,1,0.1, and 0.01) analyzed using the modified Pearson's 
statistic [definition: eq. (P)] (same input data set as for Fig. In all cases, z/ = 10^ and rrii 
was set to the true mean value of the data set. Other details as in Fig. [1|. 
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Fig. 7. — Reduced chi-square as a function of the true Poisson mean for the modified 
Pearson's statistic with the model of the zth deviate set to the true mean value of the 
parent Poisson distribution (same input data set as for Fig. Other details as in Fig. ^ 
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Fig. 8. — The cumulative distribution functions for 1000 samples of 10^ Poisson deviates 
{top to bottom: fi = 100,10,1,0.1, and 0.01) analyzed using the modified Pearson's 
statistic (same input data set as for Fig. In all cases, u = 10^ and rrii was set to the 
sample mean. Other details as in Fig. |1]. 
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Fig. 9. — Reduced chi-square as a function of the true Poisson mean for the modified 
Pearson's statistic with the model of the ith deviate set to the sample mean (same input 
data set as for Fig. Other details as in Fig. ^ 
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Fig. 10. — The cumulative distribution functions for 1000 samples of 10^ Poisson deviates 
{top to bottom: fi = 100,10,1,0.1, and 0.01) analyzed using the statistic [definition: eq. 
(^)] (same input data set as for Fig. 0). In all cases, i> = 10'^ and rrii was set to the true mean 
value of the data set. Other details as in Fig. 



10 



!f niiBij ["I'liBf i"i'niiiii|""i"i'iiiiiiij niirBf I'liiBiji 




0.001 0.01 0.1 1 10 100 
Irue mean : 



1000 



Fig. 11. — Reduced chi-square as a function of the true Poisson mean for the statistic 
statistic with the model of the ith deviate set to the true mean value of the parent Poisson 
distribution (same input data set as for Fig. ||). The solid line connecting the open squares 
in the top panel is the formula 1 + (/i — 1) [eq. 29 of PaperI]. The solid line connecting 
the open squares in the bottom panel is equation (Q). Other details as in Fig. ^. 
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Fig. 12. — The cumulative distribution functions for 1000 samples of 10^ Poisson deviates 
{top to bottom: fi = 100, 10, 1,0.1, and 0.01) analyzed using the modified statistic [defi- 
nition: eq. (0)] (same input data set as for Fig. H). In all cases, h' = 10'^ and was set to 
the true mean value of the data set. Other details as in Fig. 
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Fig. 13. — Reduced chi-square as a function of the true Poisson mean for the modified 
statistic with the model of the zth deviate set to the true mean value of the parent Poisson 
distribution (same input data set as for Fig. ^). Other details as in Fig. ^ 
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Fig. 14. — The cumulative distribution functions for 1000 samples of 10^ Poisson deviates 
{top to bottom: fi = 100, 10, 1,0.1, and 0.01) analyzed using the modified statistic (same 
input data set as for Fig. 0). In all cases, i/ = 10^ and rrii was set to the sample mean. Other 
details as in Fig. 
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Fig. 15. — Reduced chi-square as a function of the true Poisson mean for the modified 
statistic with the model of the ith deviate set to the sample mean (same input data set as 
for Fig. ^. Other details as in Fig. |^. 
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Fig. 16. — The simulated X-ray observation. A 40 photon X-ray point source is located at 
the {x,y) position of (33,26) on a background of 0.06 photons per pixel. The Point Spread 
Function is (j){x,y) = (7r/3)[l — min(r, 1)] where = (x/10)^ -|- (?//10)^. This is the same 
PSF used by Cash ( |1979| ) but with a resolution of 100 pixels per unit area. 
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Fig. 17. — Gray boxes indicate pixels with a total number of photons within a radius of 10 
pixels that are consistent (within the 99.9% upper and lower single-sided confidence limits of 
the observed photon total) with the estimated background flux level of 0.0747 photons per 
pixel. The X marks the center of the X-ray point source and the dotted circle has a radius 
of 10 pixels which is the maximum size of the PSF. Other details as in Fig. ITB. 
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Fig. 18. — Gray boxes indicate the pixels with a total number of photons within a radius of 
10 pixels that are consistent (within the 99.9% upper and lower single-sided confidence limits 
of the observed photon total) with the true background flux level of 0.06 photons per pixel. 
Dark-gray circles indicate the pixels with a total number of photons within a radius of 10 
pixels that are consistent (within the 95% upper and lower single-sided confidence limits of 
the observed photon total) with the model of a 40 photon point source at that pixel location 
on a background of 0.06 photons per pixel. Other details as in Fig. 0. 
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Fig. 19. — Dark-gray circles indicate the pixels with a total number of photons within a 
radius of 8 pixels that are consistent (within the 95% upper and lower single-sided confidence 
limits of the observed photon total) with the model of a 40 photon point source at that pixel 
location on a background of 0.06 photons per pixel. All pixels within the solid black contour 
are within the 95% confidence interval as determined by the Xju analysis of the cumulative 
radial distribution of the data within 10 pixels is compared with the cumulative radial 
distribution of a model of a 40 photon point source at that pixel location on a background 
of 0.06 photons per pixel. Note how well the 95% confidence interval of the x^m analysis of 
the cumulative radial distribution matches the region described by the circled pixels. Other 



details as in Fig. [18 
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Fig. 20. — The photon distribution of the 317 pixels within a radius of 10 pixels of 
the location (33,26) of Fig. |16| was transformed to a cumulative radial distribution (10 
1-pixel-wide bins =^ 10 degrees-of-freedom) and then compared, using the modified chi- 
square-gamma statistic, with 80 models of the observation: a 1 to 80 photon (in steps of 1 
photon) X-ray point source at (33,26) on a background of 0.06 photons per pixel (i.e., the 
true background). The 95th percentage point for the chi-square distribution with 10 degrees 
of freedom is 18.31 [i.e. P(18.31|10) = 0.95]. Assuming that x^m distributed like x^y we 
see that the upper and lower single-sided 95% confidence limits for the intensity of an X-ray 
point source at (33, 26) in Fig. |16| is 54.5 and 28.0, respectively. The true intensity of the 
X-ray source is 40 photons. 
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Fig. 21. — A data set of 10^ realizations of the the same model used to make Fig. 16 was 
created. Each sample in this data set was then analyzed using the modified chi-square-gamma 
statistic at the location (33, 26) - the true location of the simulated X-ray point source of 40 
photons on a background of 0.06 photons per pixel. All 317 pixels within a radius of 10 pixels 
(the size of the PSF) were compared to the true model value at that location and the number 
of independent degrees-of-freedom was therefore equal to the number of pixels analyzed (i.e. 
z/ = 317). Compare the cumulative distribution with the cumulative distribution function of 
a Gaussian distribution with a mean of 317 and a variance of a/2 x 317 (~ 25.2) [thick curve]. 
The numbers with errors shown on the right side give the mean and rms value for the x^m 
(top) and ideal {bottom) statistics. 




Fig. 22. — The 10^ simulations of Fig. 21 were sorted by the value of the modified chi- 
square-gamma statistic. The dark plot shows the probablity P(x^]yi|317) as a function of 
the sorted x^m values. The gray plot on the bottom shows the residuals from the ideal 
one-to-one correspondance. The predicted probabilities for the 9000th, 9500th, and 9900th 
sorted simulations were 90.2069%, 94.8666%, and 98.9497%, which agrees very well with the 
expected probabilities of 90%, 95%, and 99%, respectively. For the entire simulation, the 
mean and rms value of the resiuals is 0.0013 ±0.0038 percentage points — note that the 
residuals never exceeds 1 percent. Figures ^ and indicate that the assumption that the 
X^M statistic is distributed like the distribution was valid. 



